home *** CD-ROM | disk | FTP | other *** search
/ The AGA Experience 3 / AGA Experience Volume 3 (1997)(NFA - SAdENESS)[!].iso / software / utilities / graphics / raylab / source / intersct.c < prev    next >
C/C++ Source or Header  |  1996-09-01  |  38KB  |  1,451 lines
  1. /*
  2.     name:    intersct.c
  3.  
  4.     Intersection-routines
  5.     ---------------------
  6.  
  7.     These routines will calculate the t-parameter of a line,
  8.     which corresponds to the point where the given line intersects with
  9.     the given primitive. The routines return zero (0) if no intersection
  10.     occured.
  11.  
  12.  
  13.     This source-code is part of the RayLab 1.1 package, and it is provided
  14.     for your compiling pleasure.  You may use it, change it, re-compile it
  15.     etc., as long as nobody else but you receive the changes/compilations
  16.     that you have made!
  17.  
  18.     You may not use any part(s) of this source-code for your own productions
  19.     without the permission from the author of RayLab. Please read the legal
  20.     information found in the users documentation for RayLab for more details.
  21.  
  22. */
  23.  
  24.  
  25. #include  "defs.h"
  26. #include  "extern.h"
  27.  
  28.  
  29.  
  30.  
  31. /**********************************************************************
  32.  *
  33.  *    Try all objects for intersection, and return object number
  34.  *    and intersection information.
  35.  *
  36.  **********************************************************************/
  37.  
  38. long Intersect_All_Objects(INTERSECTION *Intersct, LINE *RayLine)
  39. {
  40.     INTERSECTION TmpIntersct;
  41.     double    t, tnew;
  42.     long    i, ObjNum;
  43.  
  44.     ObjNum=-1L;
  45.  
  46.     t=-1.0;
  47.     for(i=0;i<NumObjects;i++) {
  48.         tnew=Intersect_Object(&TmpIntersct, RayLine, ObjectArray[i]);
  49.         if((tnew>EPSILON)&&((tnew<t)||(t<=EPSILON))) {
  50.         ObjNum=i;
  51.         t=tnew;
  52.         *Intersct=TmpIntersct;
  53.         }
  54.     }
  55.  
  56.     return(ObjNum);
  57. }
  58.  
  59.  
  60.  
  61. /**********************************************************************
  62.  *
  63.  *    Get (possible) intersection with a given object
  64.  *
  65.  **********************************************************************/
  66.  
  67.  
  68. double Intersect_Object(INTERSECTION *Intrsct, LINE *Line, OBJECT *Object)
  69. {
  70.     LINE    TransformedLine;
  71.     register double    t,a,b,c;
  72.     register long    i;
  73.     VECTOR    tempv;
  74.  
  75.     if(CheckForBounding(Line, Object)==TRUE) {            /* Check for boundings, if any */
  76.       TransformedLine.Origin=Line->Origin;
  77.       TransformedLine.Direction=Line->Direction;
  78.  
  79.       if(Object->Transform.NumTransforms>0) {
  80.         for(i=Object->Transform.NumTransforms-1L;i>=0L;i--) {
  81.         switch(Object->Transform.Entry[i].Type) {
  82.             case TRANSFORM_SCALE:
  83.             tempv=Object->Transform.Entry[i].Values;
  84.             a=1.0/tempv.x; b=1.0/tempv.y; c=1.0/tempv.z;
  85.             TransformedLine.Origin.x*=a; TransformedLine.Origin.y*=b; TransformedLine.Origin.z*=c;
  86.             TransformedLine.Direction.x*=a; TransformedLine.Direction.y*=b; TransformedLine.Direction.z*=c;
  87.             break;
  88.             case TRANSFORM_MOVE:
  89.             SubVector((VECTOR *)&TransformedLine.Origin,(VECTOR *)&TransformedLine.Origin,&Object->Transform.Entry[i].Values);
  90.             break;
  91.             case TRANSFORM_ROTATE:
  92.             NegVector(&tempv,&Object->Transform.Entry[i].Values);
  93.             RevRotatePoint(&TransformedLine.Origin,&tempv);
  94.             RevRotateVector(&TransformedLine.Direction,&tempv);
  95.             break;
  96.             case TRANSFORM_NONE:
  97.             default:
  98.             break;
  99.         }
  100.         }
  101.       }
  102.  
  103.       switch(Object->ShapeType) {
  104.         case SHAPE_PLANE:
  105.         t=Intersect_Plane(Intrsct, &TransformedLine, Object);
  106.         break;
  107.         case SHAPE_SPHERE:
  108.         t=Intersect_Sphere(Intrsct, &TransformedLine, Object);
  109.         break;
  110.         case SHAPE_ELLIPSOID:
  111.         t=Intersect_Ellipsoid(Intrsct, &TransformedLine, Object);
  112.         break;
  113.         case SHAPE_TRIANGLE:
  114.         t=Intersect_Triangle(Intrsct, &TransformedLine, Object);
  115.         break;
  116.         case SHAPE_TRIANGLELIST:
  117.         t=Intersect_TriangleList(Intrsct, &TransformedLine, Object);
  118.         break;
  119.         case SHAPE_BOX:
  120.         t=Intersect_Box(Intrsct, &TransformedLine, Object);
  121.         break;
  122.         case SHAPE_DISC:
  123.         t=Intersect_Disc(Intrsct, &TransformedLine, Object);
  124.         break;
  125.         case SHAPE_CYLINDER:
  126.         t=Intersect_Cylinder(Intrsct, &TransformedLine, Object);
  127.         break;
  128.         case SHAPE_CONE:
  129.         t=Intersect_Cone(Intrsct, &TransformedLine, Object);
  130.         break;
  131.         case SHAPE_TORUS:
  132.         t=Intersect_Torus(Intrsct, &TransformedLine, Object);
  133.         break;
  134.         case SHAPE_PEAK:
  135.         t=Intersect_Peak(Intrsct, &TransformedLine, Object);
  136.         break;
  137.         case SHAPE_DEOFLASKA:
  138.         t=Intersect_Deoflaska(Intrsct, &TransformedLine, Object);
  139.         break;
  140.         case SHAPE_CSG:
  141.         t=Intersect_CSG(Intrsct, &TransformedLine, Object);
  142.         break;
  143.         default:
  144.         t=0.0;
  145.         break;
  146.       }
  147.  
  148.       if(t>EPSILON) {
  149.         if(Object->Transform.NumTransforms>0) {    /* Transform intersection to its real place/orientation in space */
  150.         for(i=0L;i<Object->Transform.NumTransforms;i++) {
  151.             tempv=Object->Transform.Entry[i].Values;
  152.             switch(Object->Transform.Entry[i].Type) {
  153.             case TRANSFORM_SCALE:
  154.                 a=tempv.x; b=tempv.y; c=tempv.z;
  155.                 Intrsct->IPoint.x*=a; Intrsct->IPoint.y*=b; Intrsct->IPoint.z*=c;
  156.                 Intrsct->Normal.x*=1.0/a; Intrsct->Normal.y*=1.0/b; Intrsct->Normal.z*=1.0/c;
  157.                 break;
  158.             case TRANSFORM_MOVE:
  159.                 Intrsct->IPoint.x+=tempv.x;
  160.                 Intrsct->IPoint.y+=tempv.y;
  161.                 Intrsct->IPoint.z+=tempv.z;
  162.                 break;
  163.             case TRANSFORM_ROTATE:
  164.                 RotatePoint(&Intrsct->IPoint,&tempv);
  165.                 RotateVector(&Intrsct->Normal,&tempv);
  166.                 break;
  167.             case TRANSFORM_NONE:
  168.             default:
  169.                 break;
  170.             }
  171.         }
  172.         }
  173.         if(Object->Invert==INVERT_TRUE) {
  174.         NegVector(&Intrsct->Normal,&Intrsct->Normal);    /* "Reverse" surface */
  175.         }
  176.       }
  177.     }
  178.  
  179.     return(t);
  180. }
  181.  
  182.  
  183.  
  184. /**********************************************************************
  185.  *
  186.  *    These are the actual intersection calculation routines
  187.  *
  188.  **********************************************************************/
  189.  
  190.  
  191. double Intersect_Plane(INTERSECTION *Intrsct, LINE *Line, OBJECT *Obj)
  192. {
  193.     register double    t,pnx,pny,pnz,a,tmp;
  194.     double    lvx,lvy,lvz,lx0,ly0,lz0;
  195.     PLANE    *Plane;
  196.  
  197.     PlaneIntersectionAttempts++;
  198.     Plane=(PLANE *)Obj->Shape;
  199.  
  200.     lvx=Line->Direction.x; lvy=Line->Direction.y; lvz=Line->Direction.z;
  201.     lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
  202.     pnx=Plane->Normal.x; pny=Plane->Normal.y; pnz=Plane->Normal.z; a=Plane->a;
  203.  
  204.     tmp=pnx*lvx+pny*lvy+pnz*lvz;
  205.     if(tmp!=0.0) {
  206.         t=-(pnx*lx0+pny*ly0+pnz*lz0-a)/tmp;
  207.         if(t>EPSILON) {
  208.         PlaneIntersections++;
  209.         Intrsct->IPoint.x=lx0+lvx*t; Intrsct->IPoint.y=ly0+lvy*t; Intrsct->IPoint.z=lz0+lvz*t;
  210.         Intrsct->Normal.x=pnx; Intrsct->Normal.y=pny; Intrsct->Normal.z=pnz;
  211.         Intrsct->Texture=&Obj->Texture;
  212.         Intrsct->Shadow=Obj->Shadow;
  213.         }
  214.         else t=0.0;
  215.     }
  216.     else t=0.0;
  217.  
  218.     return(t);
  219. }
  220.  
  221.  
  222.  
  223. double Intersect_Sphere(INTERSECTION *Intrsct, LINE *Line, OBJECT *Obj)
  224. {
  225.     double    vx,vy,vz,lx0,ly0,lz0,sx0,sy0,sz0,r;
  226.  
  227.     register double    s,sroot,vxdx,vydy,vzdz,vxsqr,vysqr,vzsqr;
  228.     double    t;
  229.     double    dx,dy,dz;
  230.     double    dxsqr,dysqr,dzsqr;
  231.     double    vsum,t1,t2;
  232.     SPHERE    *Sphere;
  233.  
  234.     SphereIntersectionAttempts++;
  235.     Sphere=(SPHERE *)Obj->Shape;
  236.  
  237.     vx=Line->Direction.x; vy=Line->Direction.y; vz=Line->Direction.z;
  238.     lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
  239.     sx0=Sphere->Centre.x; sy0=Sphere->Centre.y; sz0=Sphere->Centre.z;
  240.     r=Sphere->r;
  241.  
  242.                         /* These are used more than once, */
  243.     dx=lx0-sx0; dy=ly0-sy0; dz=lz0-sz0;    /* thus much time is saved by not */
  244.     vxdx=vx*dx; vydy=vy*dy; vzdz=vz*dz;    /* performing the same operation  */
  245.     dxsqr=dx*dx; dysqr=dy*dy; dzsqr=dz*dz;    /* twice (or even more).          */
  246.     vxsqr=vx*vx; vysqr=vy*vy; vzsqr=vz*vz;
  247.  
  248.     vsum=(vxsqr+vysqr+vzsqr);
  249.  
  250.     sroot=2*(vxdx*(vydy+vzdz)+vydy*vzdz)+r*r*vsum;
  251.     sroot-=vxsqr*(dysqr+dzsqr)+vysqr*(dxsqr+dzsqr)+vzsqr*(dxsqr+dysqr);
  252.  
  253.     t=0.0;
  254.     if(sroot>=0) {
  255.         sroot=sqrt(sroot);
  256.  
  257.         s=-(vxdx+vydy+vzdz);
  258.         t1=(s+sroot)/vsum;        /* Note: vsum!=0, always */
  259.         t2=(s-sroot)/vsum;        /* Also: t2<=t1, as sroot>=0 */
  260.  
  261.         if(t2>EPSILON) t=t2;
  262.         else if(t1>EPSILON) t=t1;
  263.         if(t>EPSILON) {
  264.             SphereIntersections++;
  265.             Intrsct->IPoint.x=lx0+vx*t; Intrsct->IPoint.y=ly0+vy*t; Intrsct->IPoint.z=lz0+vz*t;
  266.             Intrsct->Normal.x=Intrsct->IPoint.x-sx0;
  267.             Intrsct->Normal.y=Intrsct->IPoint.y-sy0;
  268.             Intrsct->Normal.z=Intrsct->IPoint.z-sz0;
  269.             Intrsct->Texture=&Obj->Texture;
  270.             Intrsct->Shadow=Obj->Shadow;
  271.         }
  272.     }
  273.  
  274.     return(t);
  275. }
  276.  
  277.  
  278.  
  279. double Intersect_Ellipsoid(INTERSECTION *Intrsct, LINE *Line, OBJECT *Obj)
  280. {
  281.     double    vx,vy,vz,lx0,ly0,lz0,ex0,ey0,ez0,a,b,c;
  282.  
  283.     register double    sroot,asqr,bsqr,csqr,vxdx,vydy,vzdz;
  284.     double    t;
  285.     double    dx,dy,dz;
  286.     double    vxsqr,vysqr,vzsqr,dxsqr,dysqr,dzsqr;
  287.     double    r,s,t1,t2;
  288.     ELLIPSOID *Ellipsoid;
  289.  
  290.     EllipsoidIntersectionAttempts++;
  291.     Ellipsoid=(ELLIPSOID *)Obj->Shape;
  292.  
  293.     vx=Line->Direction.x; vy=Line->Direction.y; vz=Line->Direction.z;
  294.     lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
  295.     ex0=Ellipsoid->Centre.x; ey0=Ellipsoid->Centre.y; ez0=Ellipsoid->Centre.z;
  296.     a=Ellipsoid->Radius.x; b=Ellipsoid->Radius.y; c=Ellipsoid->Radius.z;
  297.  
  298.     asqr=a*a; bsqr=b*b; csqr=c*c;        /* These are used more than once, */
  299.     dx=lx0-ex0; dy=ly0-ey0; dz=lz0-ez0;    /* thus much time is saved by not */
  300.     vxdx=vx*dx; vydy=vy*dy; vzdz=vz*dz;    /* performing the same operation  */
  301.     dxsqr=dx*dx; dysqr=dy*dy; dzsqr=dz*dz;    /* twice (or even more).          */
  302.     vxsqr=vx*vx; vysqr=vy*vy; vzsqr=vz*vz;
  303.  
  304.     s=-((vxdx/asqr)+(vydy/bsqr)+(vzdz/csqr));
  305.  
  306.     sroot=csqr*(vxsqr*bsqr+vysqr*asqr)+vzsqr*asqr*bsqr;
  307.     sroot+=2*(vxdx*(csqr*vydy+bsqr*vzdz)+asqr*vydy*vzdz);
  308.     sroot-=vxsqr*(csqr*dysqr+bsqr*dzsqr);
  309.     sroot-=vysqr*(csqr*dxsqr+asqr*dzsqr);
  310.     sroot-=vzsqr*(bsqr*dxsqr+asqr*dysqr);
  311.  
  312.     t=0.0;
  313.     if(sroot>=0) {
  314.         sroot=sqrt(sroot)/(a*b*c);
  315.  
  316.         r=((vxsqr/asqr)+(vysqr/bsqr)+(vzsqr/csqr));
  317.  
  318.         t1=(s+sroot)/r;
  319.         t2=(s-sroot)/r;
  320.  
  321.         if((t1>EPSILON)&&((t1<t2)||(t2<EPSILON))) {
  322.             t=t1;
  323.             EllipsoidIntersections++;
  324.         }
  325.         else if((t2>EPSILON)&&((t2<t1)||(t1<EPSILON))) {
  326.             t=t2;
  327.             EllipsoidIntersections++;
  328.         }
  329.         if(t>EPSILON) {
  330.             Intrsct->IPoint.x=lx0+vx*t; Intrsct->IPoint.y=ly0+vy*t; Intrsct->IPoint.z=lz0+vz*t;
  331.             Intrsct->Normal.x=(Intrsct->IPoint.x-ex0)/(a*a);
  332.             Intrsct->Normal.y=(Intrsct->IPoint.y-ey0)/(b*b);
  333.             Intrsct->Normal.z=(Intrsct->IPoint.z-ez0)/(c*c);
  334.             Intrsct->Texture=&Obj->Texture;
  335.             Intrsct->Shadow=Obj->Shadow;
  336.         }
  337.     }
  338.  
  339.     return(t);
  340. }
  341.  
  342.  
  343.  
  344. double Intersect_Triangle(INTERSECTION *Intrsct, LINE *Line, OBJECT *Obj)
  345. {
  346.     double    t,ttmp,pnx,pny,pnz;
  347.     double    lvx,lvy,lvz,lx0,ly0,lz0;
  348.     double    ui,vi,u[4],v[4];
  349.     register double un,vn,um,vm;
  350.     short    sign,count,n,m;
  351.     TRIANGLE *Triangle;
  352.  
  353.     TriangleIntersectionAttempts++;
  354.     Triangle=(TRIANGLE *)Obj->Shape;
  355.  
  356.     lvx=Line->Direction.x; lvy=Line->Direction.y; lvz=Line->Direction.z;
  357.     lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
  358.     pnx=Triangle->Normal.x; pny=Triangle->Normal.y; pnz=Triangle->Normal.z;
  359.  
  360.     t=0.0;
  361.     ttmp=pnx*lvx+pny*lvy+pnz*lvz;
  362.     if(ttmp!=0.0) {
  363.         ttmp=-(pnx*lx0+pny*ly0+pnz*lz0-Triangle->Offset)/ttmp;
  364.         if(ttmp>EPSILON) {
  365.         switch(Triangle->Projection) {    /* Transform to 2D coordinates */
  366.             case PROJ_YZ:
  367.             ui=ly0+ttmp*lvy;
  368.             vi=lz0+ttmp*lvz;
  369.             break;
  370.             case PROJ_XZ:
  371.             ui=lx0+ttmp*lvx;
  372.             vi=lz0+ttmp*lvz;
  373.             break;
  374.             case PROJ_XY:
  375.             ui=lx0+ttmp*lvx;
  376.             vi=ly0+ttmp*lvy;
  377.             break;
  378.         }
  379.         u[0]=Triangle->Corners[0].u-ui;
  380.         u[1]=Triangle->Corners[1].u-ui;
  381.         u[2]=Triangle->Corners[2].u-ui;
  382.         v[0]=Triangle->Corners[0].v-vi;
  383.         v[1]=Triangle->Corners[1].v-vi;
  384.         v[2]=Triangle->Corners[2].v-vi;
  385.  
  386.         u[3]=u[0]; v[3]=v[0];
  387.         count=0;            /* Determine amount of sides crossed */
  388.         for(n=0;n<3;n++) {
  389.             m=n+1;
  390.             vn=v[n]; vm=v[m];
  391.             if(vn>=0.0) sign=1;  else sign=-1;
  392.             if(vm>=0.0) sign+=1; else sign-=1;
  393.             if(sign==0) {
  394.             un=u[n]; um=u[m];
  395.             if((un>=0.0)&&(um>=0.0)) count+=1;
  396.             else if((un>=0.0)||(um>=0.0)) {
  397.                 if((um-vm*(um-un)/(vm-vn))>0.0) count+=1;
  398.             }
  399.             }
  400.         }
  401.         if((count&0x01)==1) {    /* If odd amount of sides were crossed... */
  402.             TriangleIntersections++;
  403.             t=ttmp;
  404.             Intrsct->Normal.x=pnx;
  405.             Intrsct->Normal.y=pny;
  406.             Intrsct->Normal.z=pnz;
  407.             Intrsct->IPoint.x=lx0+lvx*t;
  408.             Intrsct->IPoint.y=ly0+lvy*t;
  409.             Intrsct->IPoint.z=lz0+lvz*t;
  410.             Intrsct->Texture=&Obj->Texture;
  411.             Intrsct->Shadow=Obj->Shadow;
  412.         }
  413.         }
  414.     }
  415.  
  416.     return(t);
  417. }
  418.  
  419.  
  420.  
  421. double Intersect_TriangleList(INTERSECTION *Intrsct, LINE *Line, OBJECT *Obj)
  422. {
  423.     double    t,ttmp,pnx,pny,pnz;
  424.     double    lvx,lvy,lvz,lx0,ly0,lz0;
  425.     double    ui,vi,u[4],v[4];
  426.     register double un,vn,um,vm;
  427.     short    sign,count,n,m;
  428.     TRIANGLELIST *Triangle;
  429.  
  430.     Triangle=(TRIANGLELIST *)Obj->Shape;
  431.  
  432.     lvx=Line->Direction.x; lvy=Line->Direction.y; lvz=Line->Direction.z;
  433.     lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
  434.  
  435.     t=0.0;
  436.  
  437.     while(Triangle!=NULL) {
  438.         TriangleIntersectionAttempts++;
  439.         pnx=Triangle->Normal.x; pny=Triangle->Normal.y; pnz=Triangle->Normal.z;
  440.         ttmp=pnx*lvx+pny*lvy+pnz*lvz;
  441.         if(ttmp!=0.0) {
  442.         ttmp=-(pnx*lx0+pny*ly0+pnz*lz0-Triangle->Offset)/ttmp;
  443.         if(ttmp>EPSILON) {
  444.             switch(Triangle->Projection) {    /* Transform to 2D coordinates */
  445.             case PROJ_YZ:
  446.                 ui=ly0+ttmp*lvy;
  447.                 vi=lz0+ttmp*lvz;
  448.                 break;
  449.             case PROJ_XZ:
  450.                 ui=lx0+ttmp*lvx;
  451.                 vi=lz0+ttmp*lvz;
  452.                 break;
  453.             case PROJ_XY:
  454.                 ui=lx0+ttmp*lvx;
  455.                 vi=ly0+ttmp*lvy;
  456.                 break;
  457.             }
  458.             u[0]=Triangle->Corners[0].u-ui;
  459.             u[1]=Triangle->Corners[1].u-ui;
  460.             u[2]=Triangle->Corners[2].u-ui;
  461.             v[0]=Triangle->Corners[0].v-vi;
  462.             v[1]=Triangle->Corners[1].v-vi;
  463.             v[2]=Triangle->Corners[2].v-vi;
  464.  
  465.             u[3]=u[0]; v[3]=v[0];
  466.             count=0;            /* Determine amount of sides crossed */
  467.             for(n=0;n<3;n++) {
  468.             m=n+1;
  469.             vn=v[n]; vm=v[m];
  470.             if(vn>=0.0) sign=1;  else sign=-1;
  471.             if(vm>=0.0) sign+=1; else sign-=1;
  472.             if(sign==0) {
  473.                 un=u[n]; um=u[m];
  474.                 if((un>=0.0)&&(um>=0.0)) count+=1;
  475.                 else if((un>=0.0)||(um>=0.0)) {
  476.                 if((um-vm*(um-un)/(vm-vn))>0.0) count+=1;
  477.                 }
  478.             }
  479.             }
  480.         }
  481.  
  482.         if((count&0x01)==1) {    /* If odd amount of sides were crossed... */
  483.             if((ttmp>EPSILON)&&((t<EPSILON)||(ttmp<t))) {
  484.             TriangleIntersections++;
  485.             t=ttmp;
  486.             Intrsct->Normal.x=pnx;
  487.             Intrsct->Normal.y=pny;
  488.             Intrsct->Normal.z=pnz;
  489.             Intrsct->IPoint.x=lx0+lvx*t;
  490.             Intrsct->IPoint.y=ly0+lvy*t;
  491.             Intrsct->IPoint.z=lz0+lvz*t;
  492.             }
  493.         }
  494.         }
  495.  
  496.         Triangle=Triangle->NextTriangle;
  497.     }
  498.  
  499.     if(t>EPSILON) {
  500.         Intrsct->Texture=&Obj->Texture;
  501.         Intrsct->Shadow=Obj->Shadow;
  502.     }
  503.  
  504.     return(t);
  505. }
  506.  
  507.  
  508. double Intersect_Box(INTERSECTION *Intrsct, LINE *Line, OBJECT *Obj)
  509. {
  510.     register double    t,told,ttmp;
  511.     double    lvx,lvy,lvz,lx0,ly0,lz0;
  512.     register double    ipx,ipy,ipz;
  513.     POINT    ip;
  514.     VECTOR    norm;
  515.     BOX    *Box;
  516.  
  517.     BoxIntersectionAttempts++;
  518.     Box=(BOX *)Obj->Shape;
  519.  
  520.     told=-1.0;
  521.     lvx=Line->Direction.x; lvy=Line->Direction.y; lvz=Line->Direction.z;
  522.     lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
  523.  
  524.     if(lvz!=0.0) {
  525.         ttmp=(Box->Corners[0].z-lz0)/lvz;    /* xy-plane, z-min */
  526.         if(ttmp>EPSILON) {
  527.         ipx=lx0+lvx*ttmp; ipy=ly0+lvy*ttmp;
  528.         if((ipx>=Box->Corners[0].x)&&(ipx<=Box->Corners[1].x)&&(ipy>=Box->Corners[0].y)&&(ipy<=Box->Corners[1].y)) {
  529.             told=ttmp; ip.x=ipx; ip.y=ipy; ip.z=Box->Corners[0].z; norm.x=norm.y=0.0; norm.z=-1.0;
  530.         }
  531.         }
  532.         ttmp=(Box->Corners[1].z-lz0)/lvz;    /* xy-plane, z-max */
  533.         if((ttmp>EPSILON)&&((told<0.0)||(ttmp<told))) {
  534.         ipx=lx0+lvx*ttmp; ipy=ly0+lvy*ttmp;
  535.         if((ipx>=Box->Corners[0].x)&&(ipx<=Box->Corners[1].x)&&(ipy>=Box->Corners[0].y)&&(ipy<=Box->Corners[1].y)) {
  536.             told=ttmp; ip.x=ipx; ip.y=ipy; ip.z=Box->Corners[1].z; norm.x=norm.y=0.0; norm.z=1.0;
  537.         }
  538.         }
  539.     }
  540.     if(lvy!=0.0) {
  541.         ttmp=(Box->Corners[0].y-ly0)/lvy;    /* xz-plane, y-min */
  542.         if((ttmp>EPSILON)&&((told<0.0)||(ttmp<told))) {
  543.         ipx=lx0+lvx*ttmp; ipz=lz0+lvz*ttmp;
  544.         if((ipx>=Box->Corners[0].x)&&(ipx<=Box->Corners[1].x)&&(ipz>=Box->Corners[0].z)&&(ipz<=Box->Corners[1].z)) {
  545.             told=ttmp; ip.x=ipx; ip.z=ipz; ip.y=Box->Corners[0].y; norm.x=norm.z=0.0; norm.y=-1.0;
  546.         }
  547.         }
  548.         ttmp=(Box->Corners[1].y-ly0)/lvy;    /* xz-plane, y-max */
  549.         if((ttmp>EPSILON)&&((told<0.0)||(ttmp<told))) {
  550.         ipx=lx0+lvx*ttmp; ipz=lz0+lvz*ttmp;
  551.         if((ipx>=Box->Corners[0].x)&&(ipx<=Box->Corners[1].x)&&(ipz>=Box->Corners[0].z)&&(ipz<=Box->Corners[1].z)) {
  552.             told=ttmp; ip.x=ipx; ip.z=ipz; ip.y=Box->Corners[1].y; norm.x=norm.z=0.0; norm.y=1.0;
  553.         }
  554.         }
  555.     }
  556.     if(lvy!=0.0) {
  557.         ttmp=(Box->Corners[0].x-lx0)/lvx;    /* yz-plane, x-min */
  558.         if((ttmp>EPSILON)&&((told<0.0)||(ttmp<told))) {
  559.         ipy=ly0+lvy*ttmp; ipz=lz0+lvz*ttmp;
  560.         if((ipy>=Box->Corners[0].y)&&(ipy<=Box->Corners[1].y)&&(ipz>=Box->Corners[0].z)&&(ipz<=Box->Corners[1].z)) {
  561.             told=ttmp; ip.y=ipy; ip.z=ipz; ip.x=Box->Corners[0].x; norm.y=norm.z=0.0; norm.x=-1.0;
  562.         }
  563.         }
  564.         ttmp=(Box->Corners[1].x-lx0)/lvx;    /* yz-plane, x-max */
  565.         if((ttmp>EPSILON)&&((told<0.0)||(ttmp<told))) {
  566.         ipy=ly0+lvy*ttmp; ipz=lz0+lvz*ttmp;
  567.         if((ipy>=Box->Corners[0].y)&&(ipy<=Box->Corners[1].y)&&(ipz>=Box->Corners[0].z)&&(ipz<=Box->Corners[1].z)) {
  568.             told=ttmp; ip.y=ipy; ip.z=ipz; ip.x=Box->Corners[1].x; norm.y=norm.z=0.0; norm.x=1.0;
  569.         }
  570.         }
  571.     }
  572.  
  573.     if(told>EPSILON) {
  574.         t=told;
  575.         BoxIntersections++;
  576.         Intrsct->IPoint=ip;
  577.         Intrsct->Normal=norm;
  578.         Intrsct->Texture=&Obj->Texture;
  579.         Intrsct->Shadow=Obj->Shadow;
  580.     }
  581.     else t=0.0;
  582.  
  583.     return(t);
  584. }
  585.  
  586.  
  587. double Intersect_Disc(INTERSECTION *Intrsct, LINE *Line, OBJECT *Obj)
  588. {
  589.     register double    t;
  590.  
  591.     DiscIntersectionAttempts++;
  592.  
  593.     t=Intersect_DiscShape(Intrsct, Line, (DISC *)Obj->Shape);
  594.     if(t>EPSILON) {
  595.         Intrsct->Texture=&Obj->Texture;
  596.         Intrsct->Shadow=Obj->Shadow;
  597.     }
  598.  
  599.     return(t);
  600. }
  601.  
  602.  
  603. double Intersect_DiscShape(INTERSECTION *Intrsct, LINE *Line, DISC *Disc)
  604. {
  605.     register double    t,pnx,pny,pnz,a,tmp;
  606.     double    lvx,lvy,lvz,lx0,ly0,lz0,dx,dy,dz;
  607.     POINT    ip;
  608.  
  609.     lvx=Line->Direction.x; lvy=Line->Direction.y; lvz=Line->Direction.z;
  610.     lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
  611.     pnx=Disc->Plane.Normal.x; pny=Disc->Plane.Normal.y; pnz=Disc->Plane.Normal.z; a=Disc->Plane.a;
  612.  
  613.     t=0.0;
  614.     tmp=pnx*lvx+pny*lvy+pnz*lvz;
  615.     if(tmp!=0.0) {
  616.         tmp=-(pnx*lx0+pny*ly0+pnz*lz0-a)/tmp;
  617.         if(tmp>EPSILON) {
  618.             ip.x=lx0+lvx*tmp; ip.y=ly0+lvy*tmp; ip.z=lz0+lvz*tmp;
  619.             dx=ip.x-Disc->Centre.x;
  620.             dy=ip.y-Disc->Centre.y;
  621.             dz=ip.z-Disc->Centre.z;
  622.             if((dx*dx+dy*dy+dz*dz)<(Disc->r*Disc->r)) {
  623.                 t=tmp;
  624.                 DiscIntersections++;
  625.                 Intrsct->IPoint=ip;
  626.                 Intrsct->Normal.x=pnx; Intrsct->Normal.y=pny; Intrsct->Normal.z=pnz; 
  627.             }
  628.         }
  629.     }
  630.  
  631.     return(t);
  632. }
  633.  
  634.  
  635. double Intersect_Cylinder(INTERSECTION *Intrsct, LINE *Line, OBJECT *Obj)
  636. { 
  637.     register double lvx,lvy,l0x,l0y,lvx2plvy2;
  638.     double    rsqr,a,b,c,t,t1,t2,z1,z2;
  639.     int    sidehit;
  640.     CYLINDER *Cylinder;
  641.  
  642.     CylinderIntersectionAttempts++;
  643.     Cylinder=(CYLINDER *)Obj->Shape;
  644.  
  645.     lvx=Line->Direction.x; lvy=Line->Direction.y;
  646.     l0x=Line->Origin.x; l0y=Line->Origin.y;
  647.     rsqr=Cylinder->r; rsqr*=rsqr;
  648.     lvx2plvy2=1.0/(lvx*lvx+lvy*lvy);
  649.  
  650.     a=(rsqr-l0x*l0x-l0y*l0y)*lvx2plvy2;
  651.     b=(lvx*l0x+lvy*l0y)*lvx2plvy2;
  652.     c=a+b*b;
  653.     t=0.0;
  654.     if(c>=0.0) {
  655.         c=sqrt(c);
  656.         t1=-b-c;    /*  t1 <= t2  */
  657.         t2=-b+c;
  658.         z1=z2=Line->Origin.z;
  659.         z1+=Line->Direction.z*t1;
  660.         z2+=Line->Direction.z*t2;
  661.         a=Cylinder->h;        /* a = height of cylinder */
  662.         sidehit=FALSE;
  663.         if(t1>EPSILON) {
  664.             if((z1>=0.0)&&(z1<=a)) {
  665.             t=t1;
  666.             Intrsct->IPoint.z=z1;
  667.             sidehit=TRUE;
  668.             }
  669.         }
  670.         else if(t2>EPSILON) {
  671.             if((z2>=0.0)&&(z2<=a)) {
  672.             t=t2;
  673.             Intrsct->IPoint.z=z2;
  674.             sidehit=TRUE;
  675.             }
  676.         }
  677.         b=Line->Origin.z; c=1.0/Line->Direction.z;
  678.         t1=(0.0-b)*c;            /* t1 <=> bottom hit */
  679.         t2=(a-b)*c;            /* t2 <=> top hit */
  680.         if(t1>EPSILON) {
  681.             b=l0x+lvx*t1; c=l0y+lvy*t1;
  682.             if((b*b+c*c)>rsqr) t1=0.0;        /* Outside of the cyliunder? */
  683.         }
  684.         if(t2>EPSILON) {
  685.             if((t1>EPSILON)&&(t1<t2)) t2=0.0;
  686.             else {
  687.             b=l0x+lvx*t2; c=l0y+lvy*t2;
  688.             if((b*b+c*c)>rsqr) t2=0.0;
  689.             else t1=0.0;
  690.             }
  691.         }
  692.         if((t1>EPSILON)&&((t1<t)||(t<EPSILON))) {
  693.             t=t1;
  694.             CreatePoint(&Intrsct->IPoint,l0x+lvx*t1,l0y+lvy*t1,0.0);
  695.             CreateVector(&Intrsct->Normal,0.0,0.0,-1.0);
  696.             sidehit=FALSE;
  697.         }
  698.         else if((t2>EPSILON)&&((t2<t)||(t<EPSILON))) {
  699.             t=t2;
  700.             CreatePoint(&Intrsct->IPoint,b,c,a);
  701.             CreateVector(&Intrsct->Normal,0.0,0.0,1.0);
  702.             sidehit=FALSE;
  703.         }
  704.         if(sidehit==TRUE) {
  705.             Intrsct->Normal.x=Intrsct->IPoint.x=l0x+lvx*t;
  706.             Intrsct->Normal.y=Intrsct->IPoint.y=l0y+lvy*t;
  707.             Intrsct->Normal.z=0.0;
  708.         }
  709.         if(t>=EPSILON) {
  710.             CylinderIntersections++;
  711.             Intrsct->Texture=&Obj->Texture;
  712.             Intrsct->Shadow=Obj->Shadow;
  713.         }
  714.     }
  715.  
  716.     return(t);
  717. }
  718.  
  719.  
  720. double Intersect_Cone(INTERSECTION *Intrsct, LINE *Line, OBJECT *Obj)
  721. { 
  722.     double    lvx,lvy,lvz,l0x,l0y,l0z;
  723.     double    r1,r2,h,A,B,C,a,b,c,t,t1,t2,x,y,z;
  724.     int    sidehit;
  725.     CONE    *Cone;
  726.  
  727.     ConeIntersectionAttempts++;
  728.     Cone=(CONE *)Obj->Shape;
  729.  
  730.     t = 0.0;
  731.  
  732.     r1 = Cone->r1; r2 = Cone->r2; h = Cone->h;
  733.     lvx = Line->Direction.x; lvy = Line->Direction.y; lvz = Line->Direction.z;
  734.     l0x = Line->Origin.x; l0y = Line->Origin.y; l0z = Line->Origin.z;
  735.  
  736.     c = (r2-r1)/h;
  737.     A = (lvx*lvx + lvy*lvy - c*c*lvz*lvz);
  738.     if(A!=0.0) {
  739.         A = 1.0 / A;
  740.         B = lvx*l0x + lvy*l0y - lvz*c*(l0z*c + r1);
  741.         C = l0x*l0x + l0y*l0y - l0z*c*(l0z*c + 2.0*r1) - r1*r1;
  742.  
  743.         a = A*B;
  744.         b = a*a - A*C;
  745.         if(b>=0) {
  746.         b = sqrt(b);        /* b > 0 */
  747.         t1 = -a - b;        /* t1 < t2 */
  748.         t2 = -a + b;
  749.  
  750.         sidehit = FALSE;
  751.     /* Intersection with sides */
  752.         if(t2 > EPSILON) {
  753.             z = l0z + t2*lvz;
  754.             if((z>0.0)&&(z<h)) {
  755.             t = t2; sidehit = TRUE;
  756.             }
  757.             if(t1 > EPSILON) {
  758.             z = l0z + t1*lvz;
  759.             if((z>0.0)&&(z<h)) {
  760.                 t = t1; sidehit = TRUE;
  761.             }
  762.             }
  763.         }
  764.     /* Intersection with top/bottom */
  765.         if(lvz!=0.0) {
  766.             t1 = -l0z/lvz;  t2 = (h-l0z)/lvz;
  767.             if((t1>EPSILON)&&((t<EPSILON)||(t1<t))) {
  768.             x = l0x+lvx*t1; y = l0y+lvy*t1;
  769.             if((x*x + y*y) <= (r1*r1)) {
  770.                 t = t1; sidehit = FALSE;
  771.                 Intrsct->IPoint.x = x; Intrsct->IPoint.y = y;
  772.                 Intrsct->Normal.x = Intrsct->Normal.y = 0.0;
  773.                 Intrsct->Normal.z = -1.0;
  774.             }
  775.             }
  776.             if((t2>EPSILON)&&((t<EPSILON)||(t2<t))) {
  777.             x = l0x+lvx*t2; y = l0y+lvy*t2;
  778.             if((x*x + y*y) <= (r2*r2)) {
  779.                 t = t2; sidehit = FALSE;
  780.                 Intrsct->IPoint.x = x; Intrsct->IPoint.y = y;
  781.                 Intrsct->Normal.x = Intrsct->Normal.y = 0.0;
  782.                 Intrsct->Normal.z = 1.0;
  783.             }
  784.             }
  785.         }
  786.         if(t>EPSILON) {
  787.             ConeIntersections++;
  788.             Intrsct->IPoint.z = l0z+lvz*t;
  789.             if(sidehit == TRUE) {
  790.             x = Intrsct->IPoint.x = Intrsct->Normal.x = l0x+lvx*t;
  791.             y = Intrsct->IPoint.y = Intrsct->Normal.y = l0y+lvy*t;
  792.             Intrsct->Normal.z = -sqrt(x*x + y*y)*c;
  793.             }
  794.             Intrsct->Texture = &Obj->Texture;
  795.             Intrsct->Shadow = Obj->Shadow;
  796.         }
  797.         }
  798.     }
  799.  
  800.     return(t);
  801. }
  802.  
  803.  
  804. double Torus_Function(double t, double r0, double r1, LINE *Line)
  805. {
  806.     register double    a,b,c,f;
  807.  
  808.     a=Line->Origin.x+Line->Direction.x*t;
  809.     b=Line->Origin.y+Line->Direction.y*t;
  810.     c=Line->Origin.z+Line->Direction.z*t;
  811.  
  812.     f = sqrt(a*a + b*b) - r0;
  813.     f = f*f + c*c - r1*r1;
  814.  
  815.     return(f);
  816. }
  817.  
  818.  
  819. double Torus_Derivata(double t, double r0, LINE *Line)
  820. {
  821.     register double    a,b,c,df;
  822.  
  823.     a=Line->Origin.x+Line->Direction.x*t;
  824.     b=Line->Origin.y+Line->Direction.y*t;
  825.     c=Line->Origin.z+Line->Direction.z*t;
  826.  
  827.     df = (1.0 - r0/sqrt(a*a + b*b))*(a + b);
  828.     df *= (Line->Direction.x + Line->Direction.y);
  829.     df += c*Line->Direction.z;
  830.     df += df;
  831.  
  832.     return(df);
  833. }
  834.  
  835. double Intersect_Torus(INTERSECTION *Intrsct, LINE *Line, OBJECT *Obj)
  836. {
  837.     register double    r0,r1,tn,fn,tb,fb,diff;
  838.     double    tnp1,fnp1;
  839.     double    dx,dy,dz,x0,y0,z0;
  840.     double    cylr,tcyl1o,tcyl2o,tcyl1i,tcyl2i,tpln1,tpln2,z1,z2,pr1,pr2;
  841.     double    dx2pdy2,a,b,c,ts[6],t;
  842.     int    i,tamount,sortchange,FoundHit;
  843.     TORUS    *Torus;
  844.  
  845.     TorusIntersectionAttempts++;
  846.     Torus=(TORUS *)Obj->Shape;
  847.  
  848.     dx=Line->Direction.x; dy=Line->Direction.y; dz=Line->Direction.z;
  849.     x0=Line->Origin.x; y0=Line->Origin.y; z0=Line->Origin.z;
  850.     r0=Torus->r0; r1=Torus->r1;
  851.     dx2pdy2=1.0/(dx*dx+dy*dy);
  852.  
  853.     cylr=(r0+r1);        /* Check for outer cylinder-bounding */
  854.     a=(cylr*cylr-x0*x0-y0*y0)*dx2pdy2;
  855.     b=(dx*x0+dy*y0)*dx2pdy2;
  856.     c=a+b*b;
  857.  
  858.     t=0.0;
  859.     if(c>=0.0) {        /* If found cylinder, try the rest */
  860.         c=sqrt(c);
  861.         tcyl1o=-b-c;    /*  tcyl1o < tcyl2o */
  862.         tcyl2o=-b+c;
  863.         if(tcyl2o>EPSILON) {  /* Don't do anything if BOTH hits behind */
  864.         z1=(z0+dz*tcyl1o);
  865.         z2=(z0+dz*tcyl2o);
  866.         if(!(((z1>r1)&&(z2>r1))||((z1<-r1)&&(z2<-r1)))) {
  867.             if(fabs(z1)>r1) tcyl1o=0.0;
  868.             if(fabs(z2)>r1) tcyl2o=0.0;
  869.             cylr=(r0-r1);      /* Check for inner cylinder */
  870.             a=(cylr*cylr-x0*x0-y0*y0)*dx2pdy2;
  871.             b=(dx*x0+dy*y0)*dx2pdy2;
  872.             c=a+b*b;
  873.             if(c>=0.0) {
  874.             c=sqrt(c);
  875.             tcyl1i=-b-c;    /*  tcyl1i < tcyl2i */
  876.             tcyl2i=-b+c;
  877.             z1=(z0+dz*tcyl1i);
  878.             z2=(z0+dz*tcyl2i);
  879.             if(fabs(z1)>r1) tcyl1i=0.0;
  880.             if(fabs(z2)>r1) tcyl2i=0.0;
  881.             }
  882.             else {
  883.             tcyl1i=tcyl2i=0.0;
  884.             }
  885.             tpln1=(r1-z0)/dz;    /* Check for plane intersection */
  886.             tpln2=(-r1-z0)/dz;
  887.             a=x0+dx*tpln1; b=y0+dy*tpln1; pr1=sqrt(a*a+b*b);
  888.             a=x0+dx*tpln2; b=y0+dy*tpln2; pr2=sqrt(a*a+b*b);
  889.             if((pr1>(r0+r1))||(pr1<(r0-r1))) tpln1=0.0;
  890.             if((pr2>(r0+r1))||(pr2<(r0-r1))) tpln2=0.0;
  891.             if( (tcyl1o>EPSILON)||(tcyl2o>EPSILON)||
  892.             (tcyl1i>EPSILON)||(tcyl2i>EPSILON)||
  893.             (tpln1>EPSILON)||(tpln2>EPSILON) ) {
  894.             tamount=0;
  895.             if(fabs(tcyl1o)>EPSILON) ts[tamount++]=tcyl1o;
  896.             if(fabs(tcyl2o)>EPSILON) ts[tamount++]=tcyl2o;
  897.             if(fabs(tcyl1i)>EPSILON) ts[tamount++]=tcyl1i;
  898.             if(fabs(tcyl2i)>EPSILON) ts[tamount++]=tcyl2i;
  899.             if(fabs(tpln1)>EPSILON)  ts[tamount++]=tpln1;
  900.             if(fabs(tpln2)>EPSILON)  ts[tamount++]=tpln2;
  901.             if(((tamount&1)!=0)||(tamount>4)) {
  902.                 tamount=-1;
  903.             }
  904.             else {
  905.                 sortchange=1;
  906.                 while(sortchange>0) {  /* Sort bound-hits */
  907.                 sortchange=0;       /* (closest last)  */
  908.                 for(i=0;i<tamount-1;i++) {
  909.                     if(ts[i+1]>ts[i]) {
  910.                     a=ts[i]; ts[i]=ts[i+1]; ts[i+1]=a;
  911.                     sortchange++;
  912.                     }
  913.                 }
  914.                 }
  915.                 a=ts[1];
  916.                 if(tamount>2) {
  917.                  b=ts[2]; c=ts[3];
  918.                 }
  919.                 ts[1]=(ts[0]+a)*0.5; ts[2]=a;
  920.                 tamount++;
  921.                 if(tamount>3) {
  922.                 ts[3]=b; ts[4]=(b+c)*0.5; ts[5]=c;
  923.                 tamount++;
  924.                 }
  925.             }
  926.             FoundHit=FALSE;
  927.             while((tamount>=2)&&(FoundHit==FALSE)) {
  928.                 tb=ts[tamount-2];    /* tb is further away than... */
  929.                 tn=ts[tamount-1];    /* ...tn  */
  930.                 fb=Torus_Function(tb,r0,r1,Line);
  931.                 fn=Torus_Function(tn,r0,r1,Line);
  932.                 if(((fn<0.0)&&(fb>0.0))||((fn>0.0)&&(fb<0.0))) {
  933.                 if((tb>EPSILON)||(tn>EPSILON)) {
  934.                     i=0; diff=1.0;    /* Trapets-method */
  935.                     while((diff>NUMEPSILON)&&(i<10)) {
  936.                     tnp1=tn-fn*(tb-tn)/(fb-fn);
  937.                     diff=fabs(tnp1-tn);
  938.                     if(diff>NUMEPSILON) {
  939.                         fnp1=Torus_Function(tnp1,r0,r1,Line);
  940.                         if(((fnp1<0.0)&&(fb<0.0))||((fnp1>0.0)&&(fb>0.0))) {
  941.                         fb=fn; tb=tn;    /* Switch direction */
  942.                         }
  943.                         fn=fnp1;
  944.                         i++;
  945.                     }
  946.                     tn=tnp1;
  947.                     }
  948.                     if(tn>NUMEPSILON) {
  949.                     FoundHit=TRUE;
  950.                     }
  951.                 }
  952.                 }
  953.                 else {            /* Newton-Raphson */
  954.                 if((tamount==2)||(tamount==5)) {
  955.                     tn=tb; fn=fb;
  956.                 }
  957.                 i=0; diff=1.0;
  958.                 while((diff>NUMEPSILON)&&(i<6)) {
  959.                     diff=fn/Torus_Derivata(tn,r0,Line);
  960.                     tn-=diff;
  961.                     diff=fabs(diff);
  962.                     if(diff>NUMEPSILON) {
  963.                     fn=Torus_Function(tn,r0,r1,Line);
  964.                     i++;
  965.                     }
  966.                 }
  967.                 if((diff<=NUMEPSILON)&&(tn>NUMEPSILON)) {
  968.                     FoundHit=TRUE;
  969.                 }
  970.                 }
  971.                 if(FoundHit==FALSE) {
  972.                 tamount--;
  973.                 if(tamount==4) tamount--;
  974.                 }
  975.             }
  976.             if(FoundHit==TRUE) {
  977.                 TorusIntersections++;
  978.                 t=tn;
  979.                 Intrsct->Texture=&Obj->Texture;
  980.                 Intrsct->Shadow=Obj->Shadow;
  981.                 Intrsct->IPoint.x=x0+dx*t;
  982.                 Intrsct->IPoint.y=y0+dy*t;
  983.                 Intrsct->IPoint.z=z0+dz*t;
  984.                 a=Intrsct->IPoint.x; b=Intrsct->IPoint.y;
  985.                 c=1.0-r0/sqrt(a*a+b*b);
  986.                 Intrsct->Normal.x=a*c;
  987.                 Intrsct->Normal.y=b*c;
  988.                 Intrsct->Normal.z=Intrsct->IPoint.z;
  989.             }
  990.             }
  991.         }
  992.         }
  993.     }
  994.  
  995.     return(t);
  996. }
  997.  
  998.  
  999. double Peak_Function(double t, LINE *Line)
  1000. {
  1001.     register double x,y,z,f;
  1002.  
  1003.     x = Line->Origin.x + Line->Direction.x*t;
  1004.     y = Line->Origin.y + Line->Direction.y*t;
  1005.     z = Line->Origin.z + Line->Direction.z*t;
  1006.  
  1007.     f = x*x + y*y;
  1008.     if(z>0.0) f -= 1.0/(z*z);
  1009.     else f += 1.0/(z*z);
  1010.  
  1011.     return( f );
  1012. }
  1013.  
  1014. double Peak_Derivata(double t, LINE *Line)
  1015. {
  1016.     register double x,y,z,f;
  1017.  
  1018.     x = Line->Origin.x + Line->Direction.x*t;
  1019.     y = Line->Origin.y + Line->Direction.y*t;
  1020.     z = Line->Origin.z + Line->Direction.z*t;
  1021.  
  1022.     f = x*Line->Direction.x + y*Line->Direction.y;
  1023.     if(z>0.0) f += Line->Direction.z/(z*z*z);
  1024.     else f -= Line->Direction.z/(z*z*z);
  1025.  
  1026.     return( 2.0*f );
  1027. }
  1028.  
  1029. double Intersect_Peak(INTERSECTION *Intrsct, LINE *Line, OBJECT *Obj)
  1030. {
  1031.     double        lvx,lvy,lvz,lx0,ly0,lz0;
  1032.     double        t,t_a,t_b,f1,f2,tnp1,fnp1,fprim;
  1033.     register double    A,B,C,D,a,b;
  1034.     int        i;
  1035.  
  1036.     PeakIntersectionAttempts++;
  1037.  
  1038.     lvx = Line->Direction.x; lvy = Line->Direction.y; lvz = Line->Direction.z;
  1039.     lx0 = Line->Origin.x; ly0 = Line->Origin.y; lz0 = Line->Origin.z;
  1040.  
  1041.     a = 1.0/lvx; b = 1.0/lvy;
  1042.     t_a = -(ly0*a+lx0*b)/(lvy*a+lvx*b);    /* Closest point to z-axis */
  1043.     t_b = -lz0/lvz;                /* Intersection with z=0 */
  1044.  
  1045.  
  1046.   /* Use Newton-Raphson to find hit with "floor" */
  1047.     tnp1 = t_b; b = 1.0; i=0;
  1048.     tnp1 = 0.0;
  1049.     while((fabs(b)>NUMEPSILON)&&((i++)<6)) {
  1050.         b = Peak_Function(tnp1,Line) / Peak_Derivata(tnp1,Line);
  1051.         tnp1 -= b;
  1052.     }
  1053.     t_b = tnp1;                /* t_b = intersection with "floor" */
  1054.  
  1055.  
  1056.  /*    a = t_a;
  1057.     f1 = D + C*a;
  1058.     a *= t_a; f1 += B*a;
  1059.     a *= t_a; f1 += A*a; */            /* f1 = f(t_a) */
  1060.  
  1061.  /*    if(f1>0.0) !utanför! */
  1062.  
  1063.     t = t_b;
  1064.  
  1065.     if(t>(NUMEPSILON*4.0)) {
  1066.         PeakIntersections++;
  1067.         a = lx0 + lvx*t; b = ly0 + lvy*t;
  1068.         Intrsct->IPoint.x=a; Intrsct->Normal.x=a;
  1069.         Intrsct->IPoint.y=b; Intrsct->Normal.y=b;
  1070.         a = lz0 + lvz*t;
  1071.         Intrsct->IPoint.z=a; Intrsct->Normal.y=1.0/(a*a*a);
  1072.         NormalizeVector(&Intrsct->Normal,&Intrsct->Normal); /* Can be _HUGE_ and _SMALL_! */
  1073.         Intrsct->Texture = &Obj->Texture;
  1074.         Intrsct->Shadow = Obj->Shadow;
  1075.     }
  1076.     else t = 0.0;
  1077.  
  1078.     return(t);
  1079. }
  1080.  
  1081.  
  1082. double Intersect_Deoflaska(INTERSECTION *Intrsct, LINE *Line, OBJECT *Obj)
  1083. {
  1084.     register double    t,r0,dr,h,tn,tnp1,t1,t2,diff;
  1085.     double    dx,dy,dz,x0,y0,z0,unitdz;
  1086.     double    tcyl1,tcyl2,cylr,dx2pdy2,a,b,c,z1,z2;
  1087.     int    i, MultiplePeriods, LastCheck, FoundHit;
  1088.     DEOFLASKA    *Deo;
  1089.  
  1090.     DeoflaskaIntersectionAttempts++;
  1091.     Deo=(DEOFLASKA *)Obj->Shape;
  1092.  
  1093.     dx=Line->Direction.x; dy=Line->Direction.y; dz=Line->Direction.z;
  1094.     x0=Line->Origin.x; y0=Line->Origin.y; z0=Line->Origin.z;
  1095.     r0=Deo->r0; dr=Deo->dr; h=Deo->h;
  1096.  
  1097.     cylr=(r0+dr)*1.03;    /* Check for cylinder-bounding */
  1098.     dx2pdy2=1.0/(dx*dx+dy*dy);
  1099.  
  1100.     a=(cylr*cylr-x0*x0-y0*y0)*dx2pdy2;
  1101.     b=(dx*x0+dy*y0)*dx2pdy2;
  1102.     c=a+b*b;
  1103.  
  1104.     if(c>=0.0) {        /* If found cylinder, try numerical root-finding */
  1105.         c=sqrt(c);
  1106.         tcyl1=-b-c;        /*  tcyl1 < tcyl2 */
  1107.         tcyl2=-b+c;
  1108.         if(tcyl2>EPSILON) {    /* Don't do anything if BOTH hits behind */
  1109.         FoundHit=FALSE;
  1110.         LastCheck=FALSE;
  1111.         if(((dz*dz)/(dx*dx+dy*dy))>0.25) MultiplePeriods=TRUE;
  1112.         else MultiplePeriods=FALSE;
  1113.         t1=t2=tcyl1;
  1114.         unitdz=fabs(dz)/sqrt(dx*dx+dy*dy+dz*dz);
  1115.         while((FoundHit==FALSE)&&(LastCheck==FALSE)) {
  1116.             t1=t2;        /* Determine interval t1 <= t <= t2 containing max one root */
  1117.             if(MultiplePeriods==TRUE) {
  1118.             z1=(z0+dz*t1);
  1119.             if(dz>0.0) {
  1120.                 z2=(floor((z1+PID2)/PI)*PI)+PID2;
  1121. /*                z2=floor((z1+PI)/PI)*PI; */
  1122.                 if((z2-z1)<EPSILON) z2+=PI;
  1123.             } else {
  1124.                 z2=(floor((z1-PID2)/PI)*PI)+PID2;
  1125. /*                z2=floor(z1/PI)*PI; */
  1126.                 if((z1-z2)<EPSILON) z2-=PI;
  1127.             }
  1128.             t2=(z2-z0)/dz;
  1129.             if(t2>=tcyl2) {
  1130.                 t2=tcyl2;
  1131.                 LastCheck=TRUE;
  1132.             }
  1133.             tn=(t1+t2)*0.5;    /* Start in the middle of the interval */
  1134.             }
  1135.             else {
  1136.             if(t2<=(tcyl1+NUMEPSILON)) {
  1137.                 t1=tcyl1;
  1138.                 t2=(tcyl1+tcyl2)*0.5;
  1139.                 tn=tcyl1;
  1140.             }
  1141.             else {
  1142.                 t1=(tcyl1+tcyl2)*0.5;
  1143.                 t2=tcyl2;
  1144.                 tn=tcyl2;
  1145.                 LastCheck=TRUE;
  1146.             }
  1147.             }
  1148.  
  1149.             if(t2>EPSILON) {
  1150.             i=0; diff=1.0;
  1151.             while((diff>NUMEPSILON)&&(i<6)) {
  1152.                 a=(x0+dx*tn); b=(y0+dy*tn); c=(r0-dr*sin(z0+dz*tn));
  1153.                 diff=(a*a+b*b-c*c)/(2*a*dx+2*b*dy+2*c*dr*cos(z0+dz*tn)*dz);
  1154.                 tnp1=tn-diff;
  1155.                 diff=fabs(diff);
  1156.                 tn=tnp1;
  1157.                 i++;
  1158.             }
  1159.             if((diff<NUMEPSILON)&&(tn>(4.0*NUMEPSILON))&&(tn>=t1)&&(tn<=t2)) FoundHit=TRUE;
  1160.             else FoundHit=FALSE;
  1161.             }
  1162.         }
  1163.         if((FoundHit==TRUE)&&(tn>(4.0*NUMEPSILON))&&(tn>=tcyl1)&&(tn<=tcyl2)) {
  1164.             DeoflaskaIntersections++;
  1165.             t=tn;
  1166.             Intrsct->Texture=&Obj->Texture;
  1167.             Intrsct->Shadow=Obj->Shadow;
  1168.             Intrsct->IPoint.x=Intrsct->Normal.x=x0+dx*t;
  1169.             Intrsct->IPoint.y=Intrsct->Normal.y=y0+dy*t;
  1170.             Intrsct->IPoint.z=z0+dz*t;
  1171.             Intrsct->Normal.z=(r0-dr*sin(Intrsct->IPoint.z))*dr*cos(Intrsct->IPoint.z);
  1172.         }
  1173.         else t=0.0;
  1174.         }
  1175.         else t=0.0;
  1176.     }
  1177.     else t=0.0;
  1178.  
  1179.     return(t);
  1180. }
  1181.  
  1182.  
  1183. double Intersect_CSG(INTERSECTION *Intrsct, LINE *Line, OBJECT *Obj)
  1184. {
  1185.     register double    t, t1, t2, ttmp;
  1186.     short    operator, inside1, inside2, leftobj;
  1187.     INTERSECTION Intersct1, Intersct2;
  1188.     LINE    NewLine;
  1189.     CSG    *CsgObject;
  1190.  
  1191.     CsgIntersectionAttempts++;
  1192.     CsgObject=(CSG *)Obj->Shape;
  1193.  
  1194.     t=0.0; operator=CsgObject->Operator;
  1195.  
  1196.     switch(operator) {
  1197.  
  1198.  
  1199.     /* --- Logical AND (also used for MINUS) --- */
  1200.  
  1201.         case CSG_OP_AND:
  1202.         t=t1=t2=0.0;
  1203.         inside1=inside2=CSG_OUTSIDE;
  1204.  
  1205.         leftobj=FALSE;
  1206.         NewLine.Direction=Line->Direction;
  1207.         NewLine.Origin=Line->Origin;
  1208.         while((leftobj==FALSE)&&(inside1==CSG_OUTSIDE)) {
  1209.             ttmp=Intersect_Object(&Intersct1,&NewLine,CsgObject->Object1);
  1210.             if(ttmp>EPSILON) {
  1211.             t1+=ttmp;
  1212.             inside1=InsideObject(CsgObject->Object2,&Intersct1.IPoint);
  1213.             NewLine.Origin=Intersct1.IPoint;
  1214.             }
  1215.             else leftobj=TRUE;
  1216.         }
  1217.  
  1218.         leftobj=FALSE;
  1219.         NewLine.Origin=Line->Origin;
  1220.         while((leftobj==FALSE)&&(inside2==CSG_OUTSIDE)) {
  1221.             ttmp=Intersect_Object(&Intersct2,&NewLine,CsgObject->Object2);
  1222.             if(ttmp>EPSILON) {
  1223.             t2+=ttmp;
  1224.             inside2=InsideObject(CsgObject->Object1,&Intersct2.IPoint);
  1225.             NewLine.Origin=Intersct2.IPoint;
  1226.             }
  1227.             else leftobj=TRUE;
  1228.         }
  1229.  
  1230.         if((inside1==CSG_INSIDE)&&(inside2==CSG_OUTSIDE)) {
  1231.             t=t1;
  1232.             CopyIntersection(Intrsct,&Intersct1);
  1233.         }
  1234.         else if((inside1==CSG_OUTSIDE)&&(inside2==CSG_INSIDE)) {
  1235.             t=t2;
  1236.             CopyIntersection(Intrsct,&Intersct2);
  1237.         }
  1238.         else if((inside1==CSG_INSIDE)&&(inside2==CSG_INSIDE)) {
  1239.             if(t1<t2) {
  1240.             t=t1;
  1241.             CopyIntersection(Intrsct,&Intersct1);
  1242.             } else {
  1243.             t=t2;
  1244.             CopyIntersection(Intrsct,&Intersct2);
  1245.             }
  1246.         }
  1247.  
  1248.         break;
  1249.  
  1250.  
  1251.     /* --- Logical OR (single surface) --- */
  1252.  
  1253.         case CSG_OP_OR:
  1254.         t=t1=t2=0.0;
  1255.         inside1=inside2=CSG_INSIDE;
  1256.  
  1257.         leftobj=FALSE;
  1258.         NewLine.Direction=Line->Direction;
  1259.         NewLine.Origin=Line->Origin;
  1260.         while((leftobj==FALSE)&&(inside1==CSG_INSIDE)) {
  1261.             ttmp=Intersect_Object(&Intersct1,&NewLine,CsgObject->Object1);
  1262.             if(ttmp>EPSILON) {
  1263.             t1+=ttmp;
  1264.             inside1=InsideObject(CsgObject->Object2,&Intersct1.IPoint);
  1265.             NewLine.Origin=Intersct1.IPoint;
  1266.             }
  1267.             else leftobj=TRUE;
  1268.         }
  1269.  
  1270.         leftobj=FALSE;
  1271.         NewLine.Origin=Line->Origin;
  1272.         while((leftobj==FALSE)&&(inside2==CSG_INSIDE)) {
  1273.             ttmp=Intersect_Object(&Intersct2,&NewLine,CsgObject->Object2);
  1274.             if(ttmp>EPSILON) {
  1275.             t2+=ttmp;
  1276.             inside2=InsideObject(CsgObject->Object1,&Intersct2.IPoint);
  1277.             NewLine.Origin=Intersct2.IPoint;
  1278.             }
  1279.             else leftobj=TRUE;
  1280.         }
  1281.  
  1282.         if((inside1==CSG_OUTSIDE)&&(inside2==CSG_INSIDE)) {
  1283.             t=t1;
  1284.             CopyIntersection(Intrsct,&Intersct1);
  1285.         }
  1286.         else if((inside1==CSG_INSIDE)&&(inside2==CSG_OUTSIDE)) {
  1287.             t=t2;
  1288.             CopyIntersection(Intrsct,&Intersct2);
  1289.         }
  1290.         else if((inside1==CSG_OUTSIDE)&&(inside2==CSG_OUTSIDE)) {
  1291.             if(t1<t2) {
  1292.             t=t1;
  1293.             CopyIntersection(Intrsct,&Intersct1);
  1294.             } else {
  1295.             t=t2;
  1296.             CopyIntersection(Intrsct,&Intersct2);
  1297.             }
  1298.         }
  1299.  
  1300.         break;
  1301.  
  1302.  
  1303.     /* --- PLUS (closest hit, simple OR) --- */
  1304.  
  1305.         case CSG_OP_PLUS:
  1306.         t1=Intersect_Object(&Intersct1,Line,CsgObject->Object1);
  1307.         t2=Intersect_Object(&Intersct2,Line,CsgObject->Object2);
  1308.         if(t1>EPSILON) {
  1309.             if(t1<t2) t=t1;
  1310.             else if(t2>EPSILON) t=t2;
  1311.             else t=t1;
  1312.         }
  1313.         else if(t2>EPSILON) t=t2;
  1314.  
  1315.         if(t>EPSILON) {
  1316.             if(fabs(t-t1)<EPSILON) CopyIntersection(Intrsct,&Intersct1);
  1317.             else      CopyIntersection(Intrsct,&Intersct2);
  1318.         }
  1319.  
  1320.         break;
  1321.  
  1322.  
  1323.         default:
  1324.         break;
  1325.  
  1326.     }
  1327.  
  1328.     if(t>0) CsgIntersections++;
  1329.  
  1330.     return(t);
  1331. }
  1332.  
  1333.  
  1334.  
  1335. /**********************************************************************
  1336.  *
  1337.  *    Test for bounding intersections
  1338.  *
  1339.  **********************************************************************/
  1340.  
  1341. int CheckForBounding(LINE *Line, OBJECT *Obj)
  1342. {
  1343.     switch(Obj->BoundType) {
  1344.        case BOUND_BOX:
  1345.         return(CheckForBoundingBox(Line,(BOX *)Obj->Bound));
  1346.         break;
  1347.        case BOUND_SPHERE:
  1348.         return(CheckForBoundingSphere(Line,(SPHERE *)Obj->Bound));
  1349.         break;
  1350.        case BOUND_NONE:
  1351.        default:
  1352.         return(TRUE);
  1353.         break;
  1354.     }
  1355. }
  1356.  
  1357. int CheckForBoundingBox(LINE *Line, BOX *Box)
  1358. {
  1359.     register double    ttmp;
  1360.     double    lvx,lvy,lvz,lx0,ly0,lz0;
  1361.     register double    ipx,ipy,ipz;
  1362.  
  1363.     lvx=Line->Direction.x; lvy=Line->Direction.y; lvz=Line->Direction.z;
  1364.     lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
  1365.  
  1366.     if(lvz!=0.0) {
  1367.         ttmp=(Box->Corners[0].z-lz0)/lvz;    /* xy-plane, z-min */
  1368.         if(ttmp>EPSILON) {
  1369.         ipx=lx0+lvx*ttmp; ipy=ly0+lvy*ttmp;
  1370.         if((ipx>=Box->Corners[0].x)&&(ipx<=Box->Corners[1].x)&&(ipy>=Box->Corners[0].y)&&(ipy<=Box->Corners[1].y))
  1371.             return(TRUE);
  1372.         }
  1373.         ttmp=(Box->Corners[1].z-lz0)/lvz;    /* xy-plane, z-max */
  1374.         if(ttmp>EPSILON) {
  1375.         ipx=lx0+lvx*ttmp; ipy=ly0+lvy*ttmp;
  1376.         if((ipx>=Box->Corners[0].x)&&(ipx<=Box->Corners[1].x)&&(ipy>=Box->Corners[0].y)&&(ipy<=Box->Corners[1].y))
  1377.             return(TRUE);
  1378.         }
  1379.     }
  1380.     if(lvy!=0.0) {
  1381.         ttmp=(Box->Corners[0].y-ly0)/lvy;    /* xz-plane, y-min */
  1382.         if(ttmp>EPSILON) {
  1383.         ipx=lx0+lvx*ttmp; ipz=lz0+lvz*ttmp;
  1384.         if((ipx>=Box->Corners[0].x)&&(ipx<=Box->Corners[1].x)&&(ipz>=Box->Corners[0].z)&&(ipz<=Box->Corners[1].z))
  1385.             return(TRUE);
  1386.         }
  1387.         ttmp=(Box->Corners[1].y-ly0)/lvy;    /* xz-plane, y-max */
  1388.         if(ttmp>EPSILON) {
  1389.         ipx=lx0+lvx*ttmp; ipz=lz0+lvz*ttmp;
  1390.         if((ipx>=Box->Corners[0].x)&&(ipx<=Box->Corners[1].x)&&(ipz>=Box->Corners[0].z)&&(ipz<=Box->Corners[1].z))
  1391.             return(TRUE);
  1392.         }
  1393.     }
  1394.     if(lvx!=0.0) {
  1395.         ttmp=(Box->Corners[0].x-lx0)/lvx;    /* yz-plane, x-min */
  1396.         if(ttmp>EPSILON) {
  1397.         ipy=ly0+lvy*ttmp; ipz=lz0+lvz*ttmp;
  1398.         if((ipy>=Box->Corners[0].y)&&(ipy<=Box->Corners[1].y)&&(ipz>=Box->Corners[0].z)&&(ipz<=Box->Corners[1].z))
  1399.             return(TRUE);
  1400.         }
  1401.         ttmp=(Box->Corners[1].x-lx0)/lvx;    /* yz-plane, x-max */
  1402.         if(ttmp>EPSILON) {
  1403.         ipy=ly0+lvy*ttmp; ipz=lz0+lvz*ttmp;
  1404.         if((ipy>=Box->Corners[0].y)&&(ipy<=Box->Corners[1].y)&&(ipz>=Box->Corners[0].z)&&(ipz<=Box->Corners[1].z))
  1405.             return(TRUE);
  1406.         }
  1407.     }
  1408.  
  1409.     return(FALSE);
  1410. }
  1411.  
  1412.  
  1413.  
  1414. int CheckForBoundingSphere(LINE *Line, SPHERE *Sphere)
  1415. {
  1416.     double    vx,vy,vz,r;
  1417.     register double    sroot,vydy,vzdz,vxsqr,vysqr,vzsqr;
  1418.     double    dx,dy,dz;
  1419.     double    dxsqr,dysqr,dzsqr;
  1420.  
  1421.     vx=Line->Direction.x; vy=Line->Direction.y; vz=Line->Direction.z;
  1422.     r=Sphere->r;
  1423.  
  1424.     dx=Line->Origin.x-Sphere->Centre.x;
  1425.     dy=Line->Origin.y-Sphere->Centre.y;
  1426.     dz=Line->Origin.z-Sphere->Centre.z;
  1427.     vydy=vy*dy; vzdz=vz*dz;
  1428.     dxsqr=dx*dx; dysqr=dy*dy; dzsqr=dz*dz;
  1429.     vxsqr=vx*vx; vysqr=vy*vy; vzsqr=vz*vz;
  1430.  
  1431.     sroot=2*(vx*dx*(vydy+vzdz)+vydy*vzdz)+r*r*(vxsqr+vysqr+vzsqr);
  1432.     sroot-=vxsqr*(dysqr+dzsqr)+vysqr*(dxsqr+dzsqr)+vzsqr*(dxsqr+dysqr);
  1433.  
  1434.     if(sroot>=0.0) return(TRUE);
  1435.     else           return(FALSE);
  1436. }
  1437.  
  1438.  
  1439.  
  1440.  
  1441. /*  Misc. */
  1442.  
  1443. void CopyIntersection(INTERSECTION *i1, INTERSECTION *i2)
  1444. {
  1445.     i1->IPoint=i2->IPoint;
  1446.     i1->Normal=i2->Normal;
  1447.     i1->Texture=i2->Texture;
  1448.     i1->Shadow=i2->Shadow;
  1449. }
  1450.  
  1451.